#####State-level crude mortality analysis ######
setwd(<INSERT WD HERE>)
#State-level/national
az<-read.table("Arizona, 1979-1998.txt",header=T,as.is=T)
az<-rbind(az,read.table("Arizona, 1999-2014.txt",header=T,as.is=T))
az$state<-"AZ"
az$treat<-1
me<-read.table("Maine, 1979-1998.txt",header=T,as.is=T)
me<-rbind(me,read.table("Maine, 1999-2014.txt",header=T,as.is=T))
me$state<-"ME"
me$treat<-1
nv<-read.table("Nevada, 1979-1998.txt",header=T,as.is=T)
nv<-rbind(nv,read.table("Nevada, 1999-2014.txt",header=T,as.is=T))
nv$state<-"NV"
nv$treat<-0
nh<-read.table("New Hampshire, 1979-1998.txt",header=T,as.is=T)
nh<-rbind(nh,read.table("New Hampshire, 1999-2014.txt",header=T,as.is=T))
nh$state<-"NH"
nh$treat<-0
nm<-read.table("New Mexico, 1979-1998.txt",header=T,as.is=T)
nm<-rbind(nm,read.table("New Mexico, 1999-2014.txt",header=T,as.is=T))
nm$state<-"NM"
nm$treat<-0
ny<-read.table("New York, 1979-1998.txt",header=T,as.is=T)
ny<-rbind(ny,read.table("New York, 1999-2014.txt",header=T,as.is=T))
ny$state<-"NY"
ny$treat<-1
pa<-read.table("Pennsylvania, 1979-1998.txt",header=T,as.is=T)
pa<-rbind(pa,read.table("Pennsylvania, 1999-2014.txt",header=T,as.is=T))
pa$state<-"PA"
pa$treat<-0
us<-read.table("US, 1979-1998.txt",header=T,as.is=T)
us<-rbind(us,read.table("US, 1999-2014.txt",header=T,as.is=T))
us$state<-"US"
us$treat<--1

cmf<-rbind(az,me,nv,nh,nm,ny,pa,us)

#Clean
cmf$deaths <- as.numeric(cmf$Deaths)
cmf$pop <- as.numeric(cmf$Population)
cmf$crude.rate <- as.numeric(cmf$Crude.Rate)
cmf$age <- as.factor(cmf$Age.Group)
cmf$sex <- as.factor(cmf$Gender)
cmf$race <- as.factor(ifelse(cmf$Race.Code %in% c("1002-5","2131-1","A-PI"),"Other",cmf$Race))
cmf$year=cmf$Year
cmf$death.impute = as.integer(is.na(cmf$deaths))
cmf.slim <- subset(cmf,!is.na(pop))
cmf.slim$deaths=ifelse(cmf.slim$death.impute==1,1,cmf.slim$deaths) #Impute 1 death for suppressed

#Crude mortality - overall
overallGroup<-aggregate(cmf.slim[,c("deaths","pop")], by = list(cmf.slim$year,cmf.slim$treat), FUN = sum)
overallGroup$rate<-(overallGroup$deaths/overallGroup$pop)*100000
overallGroup

#Crude mortality - by state
overallState<-aggregate(cmf.slim[,c("deaths","pop")], by = list(cmf.slim$year,cmf.slim$state), FUN = sum)
overallState$rate<-(overallState$deaths/overallState$pop)*100000

#New York
ny<-subset(cmf.slim,state="ny")
nyAge<-aggregate(ny[,c("deaths","pop")], by=list(ny$year,ny$Age.Group), FUN=sum)
nyAge$rate<-(nyAge$deaths/nyAge$pop)*100000
nyRace<-aggregate(ny[,c("deaths","pop")], by=list(ny$year,ny$race), FUN=sum)
nyRace$rate<-(nyRace$deaths/nyRace$pop)*100000
nySR<-aggregate(ny[,c("deaths","pop")], by=list(ny$year,ny$sex,ny$race), FUN=sum)
nySR$rate<-(nySR$deaths/nySR$pop)*100000
#External causes of death
nyExt<-read.table("New York External, 1979-1998.txt",header=T,as.is=T)
nyExt<-rbind(nyExt,read.table("New York External, 1999-2014.txt",header=T,as.is=T))
nyExt$deaths <- as.numeric(nyExt$Deaths)
nyExt$pop <- as.numeric(nyExt$Population)
nyExt$crude.rate <- as.numeric(nyExt$Crude.Rate)
nyExt$age <- as.factor(nyExt$Age.Group)
nyExt$sex <- as.factor(nyExt$Gender)
nyExt$race <- as.factor(ifelse(nyExt$Race.Code %in% c("1002-5","2131-1","A-PI"),"Other",nyExt$Race))
nyExt$year=nyExt$Year
nyExt$death.impute = as.integer(is.na(nyExt$deaths))
nyExt.slim <- subset(nyExt,!is.na(pop))
nyExt.slim$deaths=ifelse(nyExt$death.impute==1,1,nyExt$deaths) #Impute 1 death for suppressed
ext<-aggregate(nyExt.slim[,c("deaths","pop")], by=list(nyExt.slim$year,nyExt.slim$sex,nyExt.slim$race), FUN=sum)
ext$rate<-(ext$deaths/ext$pop)*100000
#HIV
nyHIV<-read.table("New York HIV, 1979-1998.txt",header=T,as.is=T)
nyHIV<-rbind(nyHIV,read.table("New York HIV, 1999-2014.txt",header=T,as.is=T))
nyHIV$deaths <- as.numeric(nyHIV$Deaths)
nyHIV$pop <- as.numeric(nyHIV$Population)
nyHIV$sex <- as.factor(nyHIV$Gender)
nyHIV$race <- as.factor(ifelse(nyHIV$Race.Code %in% c("1002-5","2131-1","A-PI"),"Other",nyHIV$Race))
nyHIV$year=nyHIV$Year
nyHIV$death.impute = as.integer(is.na(nyHIV$deaths))
nyHIV.slim <- subset(nyHIV,!is.na(pop))
nyHIV.slim$deaths=ifelse(nyHIV$death.impute==1,1,nyHIV$deaths) #Impute 1 death for suppressed
hiv<-aggregate(nyHIV.slim[,c("deaths","pop")], by=list(nyHIV.slim$year,nyHIV.slim$sex,nyHIV.slim$race), FUN=sum)
hiv$rate<-(hiv$deaths/hiv$pop)*100000


#####County-level mortality analysis ######

####### Read in county-level mortality and clean ####################
setwd(<SET WD HERE>)

addState<-function(cmf, name){
  races<-c("White","Black or African American")
  prev=read.table(paste(name,", 1979-1998.txt",sep=""),header=T,as.is=T)
  prev$Race.Code<-NULL
  cur=read.table(paste(name,", 1999-2014.txt",sep=""),header=T,as.is=T)
  cur=subset(cur,Race %in% races) #Drop Asian/American Indian
  cur$Race.Code<-NULL #Drop race code to bind with Other
  curOther=read.table(paste(name," Other, 1999-2014.txt",sep=""),header=T,as.is=T) #Get pooled Other races
  curOther$Race="Other Race"
  cmf=rbind(cmf,prev,cur,curOther) #Append state to master
  return(cmf)
}

cmf = data.frame() #Initialize master
cmf=addState(cmf,"Arizona")
cmf=addState(cmf,"California")
cmf=addState(cmf,"Connecticut")
cmf=addState(cmf,"Delaware")
cmf=addState(cmf,"Illinois")
cmf=addState(cmf,"Iowa")
cmf=addState(cmf,"Kansas")
cmf=addState(cmf,"Maine")
cmf=addState(cmf,"Maryland")
cmf=addState(cmf,"Minnesota")
cmf=addState(cmf,"Missouri")
cmf=addState(cmf,"Nebraska")
cmf=addState(cmf,"Nevada")
cmf=addState(cmf,"New Hampshire")
cmf=addState(cmf,"New Jersey")
cmf=addState(cmf,"New Mexico")
cmf=addState(cmf,"New York")
cmf=addState(cmf,"Pennsylvania")
cmf=addState(cmf,"South Dakota")
cmf=addState(cmf,"Utah")
cmf=addState(cmf,"Vermont")
cmf=addState(cmf,"Virginia")
cmf=addState(cmf,"West Virginia")
cmf=addState(cmf,"Wisconsin")

cmf$deaths <- as.numeric(cmf$Deaths)
cmf$pop <- as.numeric(cmf$Population)
cmf$crude.rate <- as.numeric(cmf$Crude.Rate)
cmf$age <- as.factor(cmf$Age.Group)
cmf$sex <- as.factor(cmf$Gender)
cmf$death.impute = as.integer(is.na(cmf$deaths))
cmf<-subset(cmf,!is.na(pop)) #Remove suppressed populations

save(cmf,file="cmf_file.Rdata")
########

####Create synthetic counterfactuals and perform DD/ITS analyses#####
#load("~/Dropbox/Replication/Mortality/County-level/cmf_file.Rdata")

cmf$race <- as.factor(ifelse(cmf$Race %in% c("Other","Other Race"),"Other",cmf$Race))
cmf$year=cmf$Year
cmf.slim <- subset(cmf,!is.na(pop)&year %in% 1990:2014)
cmf.slim$deaths = ifelse(cmf.slim$death.impute==1,1,cmf.slim$deaths)
cmf.slim$st=as.integer(floor(cmf.slim$County.Code/1000))
cmf.slim$cty=as.factor(cmf.slim$County.Code-floor(cmf.slim$County.Code/1000)*1000)
cmf.slim$mort.rate = cmf.slim$deaths/cmf.slim$pop
library(plyr)
test = ddply(cmf.slim,.(County.Code, age, sex, race),
      summarise,
      count = length(mort.rate))
cmf.slim = merge(cmf.slim,test,by=c("County.Code", "age", "sex", "race"))
cmf.slim = subset(cmf.slim,count==25)
test = subset(test,count==25)
rm(test)
t.states=c(4,19,23,24,36)
c.states=c(6,9,10,17,20,27,29,31,32,33,34,35,42,46,49,50,51,54,55)
cmf.slim$treat = as.factor(ifelse(cmf.slim$st %in% t.states,1,0))
cmf.slim$tc.pairing.1 = as.factor(ifelse(cmf.slim$st %in% c(4,6,32,35,49),1,0))
cmf.slim$tc.pairing.2 = as.factor(ifelse(cmf.slim$st %in% c(19,17,20,27,29,31,46,55),1,0))
cmf.slim$tc.pairing.3 = as.factor(ifelse(cmf.slim$st %in% c(23,33,50),1,0))
cmf.slim$tc.pairing.4 = as.factor(ifelse(cmf.slim$st %in% c(24,10,51,54),1,0))
cmf.slim$tc.pairing.5 = as.factor(ifelse(cmf.slim$st %in% c(36,34,42,9),1,0))
cmf.slim$post = as.integer(ifelse(cmf.slim$tc.pairing.1==1,cmf.slim$year>2001,
                                  ifelse(cmf.slim$tc.pairing.2==1,cmf.slim$year>2005,
                                         ifelse(cmf.slim$tc.pairing.3==1,cmf.slim$year>2002,
                                                ifelse(cmf.slim$tc.pairing.4==1,cmf.slim$year>2007,
                                                       cmf.slim$year>2001)))))

cmf.slim=cmf.slim[order(cmf.slim$County.Code,cmf.slim$age,cmf.slim$sex,cmf.slim$race,cmf.slim$year),]
treat.ctys=data.frame()
library(parcor)
library(doMC)
registerDoMC(cores=4)

## ARIZONA
for(i in 1:5){
  for(j in 1:2){
    for(k in 1:3){
      TREAT=subset(cmf.slim,treat==1&tc.pairing.1==1
               &age==levels(cmf.slim$age)[i]
               &sex==levels(cmf.slim$sex)[j]
               &race==levels(cmf.slim$race)[k])
      CONTROL=subset(cmf.slim,treat==0&tc.pairing.1==1
                 &age==levels(age)[i]
                 &sex==levels(sex)[j]
                 &race==levels(race)[k])
      TREAT.mat = as.matrix(subset(TREAT,post==0)[,"mort.rate"])
      wgt = as.matrix(subset(TREAT,post==0)[,"pop"])
      t.len = length(unique(subset(TREAT,post==0)[,"year"]))
      n.ctys = length(TREAT.mat)/t.len
      m=as.matrix(subset(reshape(CONTROL[CONTROL$post==0,c("year","st","cty","mort.rate")],timevar="year",
                                 idvar=c("st","cty"),direction="wide"),select=-c(st,cty)))
      CONTROL.mat=matrix(rep(m,n.ctys),ncol=nrow(m),byrow=T)
      m=as.matrix(subset(reshape(CONTROL[,c("year","st","cty","mort.rate")],timevar="year",
                                 idvar=c("st","cty"),direction="wide"),select=-c(st,cty)))
      CONTROL.all.mat=matrix(rep(m,n.ctys),ncol=nrow(m),byrow=T)
      print(paste(i,j,k))
      W = cv.glmnet(x=CONTROL.mat,y=TREAT.mat,weights=wgt,parallel=T)
      treat.ctys=rbind(treat.ctys,cbind(TREAT,synthetic.control=as.vector(predict(W,CONTROL.all.mat,s=W$lambda.min))))
    }
  }
}
## IOWA
for(i in 1:5){
  for(j in 1:2){
    for(k in 1:3){
      TREAT=subset(cmf.slim,treat==1&tc.pairing.2==1
                   &age==levels(cmf.slim$age)[i]
                   &sex==levels(cmf.slim$sex)[j]
                   &race==levels(cmf.slim$race)[k])
      CONTROL=subset(cmf.slim,treat==0&tc.pairing.2==1
                     &age==levels(age)[i]
                     &sex==levels(sex)[j]
                     &race==levels(race)[k])
      
      TREAT.mat = as.matrix(subset(TREAT,post==0)[,"mort.rate"])
      wgt = as.matrix(subset(TREAT,post==0)[,"pop"])
      t.len = length(unique(subset(TREAT,post==0)[,"year"]))
      n.ctys = length(TREAT.mat)/t.len
      m=as.matrix(subset(reshape(CONTROL[CONTROL$post==0,c("year","st","cty","mort.rate")],timevar="year",
                                 idvar=c("st","cty"),direction="wide"),select=-c(st,cty)))
      CONTROL.mat=matrix(rep(m,n.ctys),ncol=nrow(m),byrow=T)
      m=as.matrix(subset(reshape(CONTROL[,c("year","st","cty","mort.rate")],timevar="year",
                                 idvar=c("st","cty"),direction="wide"),select=-c(st,cty)))
      CONTROL.all.mat=matrix(rep(m,n.ctys),ncol=nrow(m),byrow=T)
      print(paste(i,j,k))
      W = cv.glmnet(x=CONTROL.mat,y=TREAT.mat,weights=wgt,parallel=T)
      treat.ctys=rbind(treat.ctys,cbind(TREAT,synthetic.control=as.vector(predict(W,CONTROL.all.mat,s=W$lambda.min))))    }
  }
}
## MAINE
for(i in 1:5){
  for(j in 1:2){
    for(k in 1:3){
      TREAT=subset(cmf.slim,treat==1&tc.pairing.3==1
                   &age==levels(cmf.slim$age)[i]
                   &sex==levels(cmf.slim$sex)[j]
                   &race==levels(cmf.slim$race)[k])
      CONTROL=subset(cmf.slim,treat==0&tc.pairing.3==1
                     &age==levels(age)[i]
                     &sex==levels(sex)[j]
                     &race==levels(race)[k])
      
      TREAT.mat = as.matrix(subset(TREAT,post==0)[,"mort.rate"])
      wgt = as.matrix(subset(TREAT,post==0)[,"pop"])
      t.len = length(unique(subset(TREAT,post==0)[,"year"]))
      n.ctys = length(TREAT.mat)/t.len
      m=as.matrix(subset(reshape(CONTROL[CONTROL$post==0,c("year","st","cty","mort.rate")],timevar="year",
                                 idvar=c("st","cty"),direction="wide"),select=-c(st,cty)))
      CONTROL.mat=matrix(rep(m,n.ctys),ncol=nrow(m),byrow=T)
      m=as.matrix(subset(reshape(CONTROL[,c("year","st","cty","mort.rate")],timevar="year",
                                 idvar=c("st","cty"),direction="wide"),select=-c(st,cty)))
      CONTROL.all.mat=matrix(rep(m,n.ctys),ncol=nrow(m),byrow=T)
      print(paste(i,j,k))
      W = cv.glmnet(x=CONTROL.mat,y=TREAT.mat,weights=wgt,parallel=T)
      treat.ctys=rbind(treat.ctys,cbind(TREAT,synthetic.control=as.vector(predict(W,CONTROL.all.mat,s=W$lambda.min))))    }
  }
}
## MARYLAND
for(i in 1:5){
  for(j in 1:2){
    for(k in 1:3){
      TREAT=subset(cmf.slim,treat==1&tc.pairing.4==1
                   &age==levels(cmf.slim$age)[i]
                   &sex==levels(cmf.slim$sex)[j]
                   &race==levels(cmf.slim$race)[k])
      CONTROL=subset(cmf.slim,treat==0&tc.pairing.4==1
                     &age==levels(age)[i]
                     &sex==levels(sex)[j]
                     &race==levels(race)[k])
      
      TREAT.mat = as.matrix(subset(TREAT,post==0)[,"mort.rate"])
      wgt = as.matrix(subset(TREAT,post==0)[,"pop"])
      t.len = length(unique(subset(TREAT,post==0)[,"year"]))
      n.ctys = length(TREAT.mat)/t.len
      m=as.matrix(subset(reshape(CONTROL[CONTROL$post==0,c("year","st","cty","mort.rate")],timevar="year",
                                 idvar=c("st","cty"),direction="wide"),select=-c(st,cty)))
      CONTROL.mat=matrix(rep(m,n.ctys),ncol=nrow(m),byrow=T)
      m=as.matrix(subset(reshape(CONTROL[,c("year","st","cty","mort.rate")],timevar="year",
                                 idvar=c("st","cty"),direction="wide"),select=-c(st,cty)))
      CONTROL.all.mat=matrix(rep(m,n.ctys),ncol=nrow(m),byrow=T)
      print(paste(i,j,k))
      W = cv.glmnet(x=CONTROL.mat,y=TREAT.mat,weights=wgt,parallel=T)
      treat.ctys=rbind(treat.ctys,cbind(TREAT,synthetic.control=as.vector(predict(W,CONTROL.all.mat,s=W$lambda.min))))    }
  }
}
## NEW YORK
for(i in 1:5){
  for(j in 1:2){
    for(k in 1:3){
      TREAT=subset(cmf.slim,treat==1&tc.pairing.5==1
                   &age==levels(cmf.slim$age)[i]
                   &sex==levels(cmf.slim$sex)[j]
                   &race==levels(cmf.slim$race)[k])
      CONTROL=subset(cmf.slim,treat==0&tc.pairing.5==1
                     &age==levels(age)[i]
                     &sex==levels(sex)[j]
                     &race==levels(race)[k])
      
      TREAT.mat = as.matrix(subset(TREAT,post==0)[,"mort.rate"])
      wgt = as.matrix(subset(TREAT,post==0)[,"pop"])
      t.len = length(unique(subset(TREAT,post==0)[,"year"]))
      n.ctys = length(TREAT.mat)/t.len
      m=as.matrix(subset(reshape(CONTROL[CONTROL$post==0,c("year","st","cty","mort.rate")],timevar="year",
                                 idvar=c("st","cty"),direction="wide"),select=-c(st,cty)))
      CONTROL.mat=matrix(rep(m,n.ctys),ncol=nrow(m),byrow=T)
      m=as.matrix(subset(reshape(CONTROL[,c("year","st","cty","mort.rate")],timevar="year",
                                 idvar=c("st","cty"),direction="wide"),select=-c(st,cty)))
      CONTROL.all.mat=matrix(rep(m,n.ctys),ncol=nrow(m),byrow=T)
      print(paste(i,j,k))
      W = cv.glmnet(x=CONTROL.mat,y=TREAT.mat,weights=wgt,parallel=T)
      treat.ctys=rbind(treat.ctys,cbind(TREAT,synthetic.control=as.vector(predict(W,CONTROL.all.mat,s=W$lambda.min))))    }
  }
}

synthetic.dataset.treat = subset(treat.ctys,select=c("st","cty","age","sex","race","tc.pairing.1",
                                                     "tc.pairing.2","tc.pairing.3","tc.pairing.4",
                                                     "tc.pairing.5","year","pop","mort.rate","post"))
synthetic.dataset.control = subset(treat.ctys,select=c("st","cty","age","sex","race","tc.pairing.1",
                                                     "tc.pairing.2","tc.pairing.3","tc.pairing.4",
                                                     "tc.pairing.5","year","pop","synthetic.control","post"))
names(synthetic.dataset.control)[13]="mort.rate"
synthetic.dataset.treat$treat=1
synthetic.dataset.control$treat=0
synthetic.dataset = rbind(synthetic.dataset.treat,synthetic.dataset.control)
synthetic.dataset$control = 1-synthetic.dataset$treat
AZ <- lm(mort.rate ~ treat:as.factor(year) + control:as.factor(year) - 1,
         data=synthetic.dataset[synthetic.dataset$tc.pairing.1==1,],
         weights=pop)
summary(AZ)
IA <- lm(mort.rate ~ treat:as.factor(year) + control:as.factor(year) - 1,
         data=synthetic.dataset[synthetic.dataset$tc.pairing.2==1,],
         weights=pop)
summary(IA)
ME <- lm(mort.rate ~ treat:as.factor(year) + control:as.factor(year) - 1,
         data=synthetic.dataset[synthetic.dataset$tc.pairing.3==1,],
         weights=pop)
summary(ME)
MD <- lm(mort.rate ~ treat:as.factor(year) + control:as.factor(year) - 1,
         data=synthetic.dataset[synthetic.dataset$tc.pairing.4==1,],
         weights=pop)
summary(MD)
NY <- lm(mort.rate ~ treat:as.factor(year) + control:as.factor(year) - 1,
         data=synthetic.dataset[synthetic.dataset$tc.pairing.5==1,],
         weights=pop)
summary(NY)

synthetic.dataset$time = as.integer(ifelse(synthetic.dataset$tc.pairing.1==1,synthetic.dataset$year-2001,
                                  ifelse(synthetic.dataset$tc.pairing.2==1,synthetic.dataset$year-2005,
                                         ifelse(synthetic.dataset$tc.pairing.3==1,synthetic.dataset$year-2002,
                                                ifelse(synthetic.dataset$tc.pairing.4==1,synthetic.dataset$year-2007,
                                                       synthetic.dataset$year-2001)))))
save(synthetic.dataset,file="~/Dropbox/Replication/AnalyticDatasets/synthetic_data.Rdata")

pool3.data <- subset(synthetic.dataset,(tc.pairing.1==1|tc.pairing.3==1|tc.pairing.5==1)&abs(time)<=11)
pool3.data$time2 <- pool3.data$time^2

pooled3 <- lm(mort.rate ~ treat*post, data=pool3.data, weights=pop)
summary(pooled3)
az.model <- lm(mort.rate ~ treat*post, data=pool3.data[pool3.data$tc.pairing.1==1,], weights=pop)
summary(az.model)
me.model <- lm(mort.rate ~ treat*post, data=pool3.data[pool3.data$tc.pairing.3==1,], weights=pop)
summary(me.model)
ny.model <- lm(mort.rate ~ treat*post, data=pool3.data[pool3.data$tc.pairing.5==1,], weights=pop)
summary(ny.model)

pooled3.its <- lm(mort.rate ~ time + post:time + post:treat:time,
                  data=pool3.data,
                  weights=pop)
summary(pooled3.its)
az.model.its <- lm(mort.rate ~ time + post:time + post:treat:time,
                   data=pool3.data[pool3.data$tc.pairing.1==1,], weights=pop)
summary(az.model.its)
me.model.its <- lm(mort.rate ~ time + post:time + post:treat:time,
                   data=pool3.data[pool3.data$tc.pairing.3==1,], weights=pop)
summary(me.model.its)
ny.model.its <- lm(mort.rate ~ time + post:time + post:treat:time,
                   data=pool3.data[pool3.data$tc.pairing.5==1,], weights=pop)
summary(ny.model.its)

pooled3.its2 <- lm(mort.rate ~ treat:time + treat:time2 + post:treat:time + post:treat:time2 +
                control:time + control:time2,
              data=pool3.data,
              weights=pop)
summary(pooled3.its2)
az.model.its2 <- lm(mort.rate ~ treat:time + treat:time2 + post:treat:time + post:treat:time2 +
                    control:time + control:time2,
                  data=pool3.data[pool3.data$tc.pairing.1==1,], weights=pop)
summary(az.model.its2)
me.model.its2 <- lm(mort.rate ~ treat:time + treat:time2 + post:treat:time + post:treat:time2 +
                     control:time + control:time2,
                   data=pool3.data[pool3.data$tc.pairing.3==1,], weights=pop)
summary(me.model.its2)
ny.model.its2 <- lm(mort.rate ~ treat:time + treat:time2 + post:treat:time + post:treat:time2 +
                     control:time + control:time2,
                   data=pool3.data[pool3.data$tc.pairing.5==1,], weights=pop)
summary(ny.model.its2)

pool5.data <- subset(synthetic.dataset,time>=-11&time<=7)
pool5.data$time2 <- pool5.data$time^2


pooled5.dd <- lm(mort.rate ~ treat*post, data=pool5.data, weights=pop)
summary(pooled5.dd)
ia.dd <- lm(mort.rate ~ treat*post, data=pool5.data[pool5.data$tc.pairing.2==1,], weights=pop)
summary(ia.dd)
md.dd <- lm(mort.rate ~ treat*post, data=pool5.data[pool5.data$tc.pairing.4==1,], weights=pop)
summary(md.dd)

pooled5.its <- lm(mort.rate ~ time + post:time + post:treat:time, data=pool5.data, weights=pop)
summary(pooled5.its)
ia.its <- lm(mort.rate ~ time + post:time + post:treat:time,
             data=pool5.data[pool5.data$tc.pairing.2==1,], weights=pop)
summary(ia.its)
md.its <- lm(mort.rate ~ time + post:time + post:treat:time,
             data=pool5.data[pool5.data$tc.pairing.4==1,], weights=pop)
summary(md.its)